#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<math.h>
#include<ctype.h>
#include <mpi.h>
#include <omp.h>


/* Inicializa uma matriz */
double **initMatriz(int ordem){
	int i;
  
	/* Cria e aloca espaço para a matriz */
	double **matrix;
	matrix = malloc(ordem*sizeof(double *));
	for (i = 0; i < ordem; i++)
		matrix[i] = malloc(ordem*sizeof(double)); 
  
	return matrix;
}

/* Inicializa um vetor */
double *initVector(int ordem){
	/* Cria e aloca espaço para o vetor */
	double *vector;
	vector = malloc(ordem*sizeof(double));
  
	return vector;
}

/* Libera os espaços em memoria alocados pelas matrizes */
void free_matrix(double **mat, int rows){
	int i;
	for(i=0; i<rows; i++){
		free(mat[i]);
	}
	free(mat);
}

/* Multiplica uma Matriz por um Vetor, apenas processando as linhas da matriz que vão desde posInit até posEnd */
void multiMat_Vec(double *result, double **mat, double *vec, int ordem, int posInit, int posEnd){
	int i, j;
	double soma;
  
	for(i = posInit; i < posEnd; i++){
		soma = 0;
		for(j = 0; j < ordem; j++){
			soma += mat[i][j] * vec[j];
		} 
		result[i - posInit] = soma;
	}  
}

/* Verifica o erro */
int verificaErro(double *X, double **MA, double *MB, int ordem, int linha_teste, double erro){
	int i;
	double result = 0;
  
	for(i = 0; i < ordem; i++){
		result += MA[linha_teste][i] * X[i];
	}
 
	if(fabs(result - MB[linha_teste]) < erro){ 
		return 0;
	}    
	return 1;
}


int main(int argc, char *argv[]){
	int i, j, iteracao;
	double result;

	/* Vetor para passar para a funcao gather quantos elementos cada processo espera e o deslocamento*/
	int *nro_elementos;
	int *deslocamento;

	/* Variáveis para registrar o tempo de início e término de execução */
	double start, end;
  
	/* Ordem da matriz */
	int J_ORDER;
  
	/* Fila a ser avaliada */
	int J_ROW_TEST;
  
	/* Erro permitido */
	double J_ERROR;
  
	/* Número máximo de iterações */
	int J_ITE_MAX;
  
	/* Matrizes utilizadas */
	double **MA, *MB, **LR, *b;
  
	/* Vetor de resposta */
	double *X, *Xaux, *X_calculado;
	
	/* Ponteiro para o arquivo de entrada */
	FILE *fp;
	
	/* Offset para posicoes do vetor */
	int posInit, posEnd;

	/* Variáveis do MPI */
	int rank, nthreads;
	
	/* Variável para verificação da inicialização do MPI e das leituras do arquivo (EOF) */
	int ret = 0;

	/* Inicializa o MPI */
	ret = MPI_Init(&argc, &argv);
	if(ret != MPI_SUCCESS){
		printf("Erro na inicialização da execução paralela com o MPI");
		return -1;
	}

	/* Funções para pegar o número de threads e o rank do processo */
	MPI_Comm_rank(MPI_COMM_WORLD, &rank);
	MPI_Comm_size(MPI_COMM_WORLD, &nthreads);
	
	/* Registra o tempo de início */
	start = MPI_Wtime();

  
	/* Abre o arquivo, de acordo com o parâmetro de entrada */
	if(argc <= 1){
		printf("Numero insuficiente de parametros. Informe o nome do arquivo.\n");
		return -1;
	}

	fp = fopen(argv[1], "r");
	if(fp < 0){
		printf("Arquivo não encontrado: \"%s\"\n",argv[1]);
		return -1;
	}

	/* Processa paralelamente a leitura do arquivo (carregando MA, MB, ...) e a criação e preparação de variaveis auxiliares */
	#pragma omp sections
	{
		/* Lê os 4 primeiros parâmetros */
		fscanf(fp, "%d%d%lf%d", &J_ORDER, &J_ROW_TEST, &J_ERROR, &J_ITE_MAX);
		
		#pragma omp section
		{		 
			/* Cria a matriz A, com o parâmetro fornecido */
			MA = initMatriz(J_ORDER);
		  
			/* Lê matriz MA (ainda não otimizado) */
			for( i = 0; i < J_ORDER; i++ ){
				for( j = 0; j < J_ORDER; j++ ){
					fscanf(fp, "%lf", &(MA[i][j]) );
				}
			}

			/* Cria o vetor B, com o parâmetro fornecido */
			MB = initVector(J_ORDER);

			/* Lê vetor MB (ainda não otimizado) */
			for( i = 0; i < J_ORDER; i++ ){
				fscanf(fp, "%lf", &(MB[i]) );
			}

			/* Fecha o arquivo */
			fclose(fp);

			/* Passos para a matriz LR, que é a soma da matriz L com a matriz R, invertendo o sinal(simplificação do algoritmo explicada no relatório) */
		  	LR = initMatriz(J_ORDER);
		  	#pragma omp parallel for default(shared) private(i,j)
			for( i = 0; i < J_ORDER; i++ ){
				for( j = 0; j < J_ORDER; j++ ){
					if(i > j) LR[i][j] =  - MA[i][j]/MA[i][i]; /* Matriz L */
					else if(i < j) LR[i][j] = - MA[i][j]/MA[i][i]; /* Matriz R */
					else LR[i][j] = 0.0; /* Diagonal zero */
				}
			}
		  
			/* Cria o vetor b, que é são os elementos de MD divididos pelos elementos da diagonal de MA */
			b = initVector(J_ORDER);
			#pragma omp parallel for default(shared) private(i)
			for( i = 0; i < J_ORDER; i++ ){
				b[i] = MB[i]/MA[i][i];
			}
		}

		#pragma omp section
		{
			/* Prepara o vetor x */
			X = initVector(J_ORDER);
			
			/* Inicializa com valores iniciais. No caso, coloquei zeros */
			#pragma omp parallel for default(shared) private(i)
			for( i = 0; i < J_ORDER; i++ ){
				X[i] = 0.0;
			}

			/* Seta os intervalos de processamento das linhas da matriz por processo */
			posInit = (rank) * (int)floor(J_ORDER/nthreads);
			posEnd = (rank+1) * (int)floor(J_ORDER/nthreads);

			/* Se for o último processo, seu posEnd deve ser a última linha da matriz */
			if(rank == nthreads - 1){
			  posEnd = J_ORDER;
			}

			X_calculado = initVector(posEnd - posInit);
			
			/* Inicializa com valores iniciais. No caso, coloquei zeros */
			#pragma omp parallel for default(shared) private(i)
			for( i = 0; i < posEnd - posInit; i++ ){
				X_calculado[i] = 0.0;
			}
		  
			/* Vetor auxiliar para computar o resultado intermediário de LR*X */
			Xaux = initVector(posEnd - posInit);

			/* Arruma parâmetros necessários para a função MPI_Allgatherv */
			nro_elementos = (int*) malloc (sizeof(int)*nthreads);
			deslocamento = (int*) malloc (sizeof(int)*nthreads);
			for(i=0;i<nthreads;i++){ 
				nro_elementos[i] = (int) floor(J_ORDER/nthreads); 
				deslocamento[i] = i*nro_elementos[i]; 
			}
			nro_elementos[nthreads-1] = J_ORDER - ((int)floor(J_ORDER/nthreads))*(nthreads-1);
		}
	}
		
	/* Começa a iteração do método */
	for(iteracao = 0; iteracao < J_ITE_MAX ; iteracao++){

		/* Faz a operação -(L + R)*X{k}. Neste caso, cada processo trabalhará em uma faixa de linhas da matriz */
		multiMat_Vec(Xaux, LR, X, J_ORDER, posInit, posEnd);

		/* Calcula o X{k+1} */
		for(i = 0; i + posInit < posEnd; i++){
			X_calculado[i] = Xaux[i] + b[i+posInit];
		}

		/* Junta as partes calculada em cada processador para formar o Xn-1 */
		MPI_Allgatherv(X_calculado,posEnd-posInit,MPI_DOUBLE,X,nro_elementos, deslocamento,MPI_DOUBLE,MPI_COMM_WORLD);
		
		/*Verifica se ja chegou no esperado*/
		if (verificaErro(X, MA, MB, J_ORDER, J_ROW_TEST, J_ERROR) == 0){
			break;
		}
	}

	/* Registra o tempo de término */
	end = MPI_Wtime();
	
	/* Espera todas as threads terminarem */
	MPI_Barrier(MPI_COMM_WORLD);

	/* Finaliza a região paralela */
	MPI_Finalize();


	if(rank == 0){
		/* Imprime os resultados */
		printf("----------------------------------------\n");
		printf("Iterations = %d\n", iteracao);
		
		result = 0;
		for(i = 0; i < J_ORDER; i++){
			result += MA[J_ROW_TEST][i] * X[i];
		}
		printf("RowTest: 2 => [%lf] =? %lf\n", result, MB[J_ROW_TEST]);  
		printf("Time = %f seconds\n", end-start);
		printf("----------------------------------------\n");
		
		/*=(Adicional para ver os resultados)=======================================*/
		/* Registra a execução em um arquivo de estatísticas */
		fp = fopen("statistics.txt", "a");
		if(fp < 0){
		  printf("Nao foi possivel adicionar resultado ao arquivo de estatisticas.\n");
		  return -1;
		}
		fprintf(fp, "----------------------------------------\n");
		fprintf(fp, "Method: Parallel. Number of nodes: %d\n", nthreads);
		fprintf(fp, "Matriz order: %d\n", J_ORDER);
		fprintf(fp, "Iterations = %d\n", iteracao);
		fprintf(fp, "RowTest: 2 => [%lf] =? %lf\n", result, MB[J_ROW_TEST]);  
		fprintf(fp, "Time = %f seconds\n", end-start);
		fclose(fp);
		
		/* Agora registra o vetor X em um arquivo específico (este não será incremental) */
		fp = fopen("resultX.txt", "w");
		if(fp < 0){
		  printf("Nao foi possivel criar um arquivo para armazenar o vetor X.\n");
		  return -1;
		}
		fprintf(fp, "vectorX = {\n");
		for( i = 0; i < J_ORDER; i++ ){
		  fprintf(fp, "      %lf,\n", X[i]);
		}
		fprintf(fp, "}\n");
		fclose(fp);
		/*======================================================================*/
	}

 	/* Libera da memória as matrizes e vetores usados */
	free_matrix(MA, J_ORDER);
	free_matrix(LR, J_ORDER);
	free(X_calculado);
	free(MB);
	free(b);
	free(X);
	free(Xaux);
	  	
	return 0;
}
